---
title: "Substance P Elisa/ MDA"
author: "Adam Bindas"
date: '2022-05-09'
output:
  html_document:
   code_folding: hide
---

```{r setup, include=FALSE}
# Load packages
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
library(gt)
library(tidyverse)
library(ggplot2)
library(pscl)
library(boot)
library(moments)
library(lmtest)
library(dplyr) 
library(ggpubr)
library(GLMMadaptive)
library(MASS)
library(RColorBrewer)
library(bbmle)
library(emmeans)
library(lmerTest)
library(lme4)

# Import Files
subp0 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Sub P Elisa/Substance P Elisa Results.csv", header=TRUE, stringsAsFactors=FALSE)
```

```{r inspection}
#Define Variables
subp0$Groups.f <- factor(subp0$group, levels = c(
  "Control", "But", "PFF",  "PFF+But", "LPS","PFF+LPS","4XPFF","4XPFF+But"))
subp0$m.id <- factor(subp0$m.id)
subp0$rat.no <- factor(subp0$rat)
subp0$plate <- factor(subp0$plate)
subp0$rerun <- factor(subp0$rerun)
subp0$qual.rem <- factor(subp0$qual.rem)
subp0$m.id <- factor(subp0$m.id)
subp0$unique.id <- factor(subp0$unique.id)

subp <- subset(subp0, qual.rem == 0)
```

```{r plots}
# Sample Size Calculation
n.count.all <- aggregate(data = subp,
                          rat ~ Groups.f,
                          function(rat) length(unique(rat)))
n.count.all

m.count.all <- aggregate(data = subp,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
m.count.all

# Combine values to a table
subp.replicates <- cbind(n.count.all[,2],m.count.all[,2])
colnames(subp.replicates) <- c("N","m")
rownames(subp.replicates) <-m.count.all[,1]
subp.replicates

# Plot Values
plot.subp1 <- ggplot(subp, aes(x=Groups.f, y=sub.p, fill=Groups.f)) +
  ylab("Substance P [pg/mL]") + 
  xlab("Groups") + 
  labs(
    title = "Substance P levels across culture groups",
    subtitle = "N = 4-7, 2 Plates")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = subp, aes(x = Groups.f, y = sub.p),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot.subp1

plot.subp2 <- ggplot(subp, aes(x=Groups.f, y=sub.p, fill=rat.no)) +
  ylab("Substance P [pg/mL]") + 
  xlab("Groups") + 
  labs(
    title = "Substance P levels across culture groups",
    subtitle = "N = 4-7, 2 Plates")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = subp, aes(x = Groups.f, y = sub.p,
                                groups = interaction(Groups.f, rat.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1)

plot.subp2

subp$low <- subp$sub.p < 20
subp1 <- subset(subp, low == 0)

m.count.all2 <- aggregate(data = subp1,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
m.count.all2

plot.subp3 <- ggplot(subp1, aes(x=Groups.f, y=sub.p, fill=Groups.f)) +
  ylab("Substance P [pg/mL]") + 
  xlab("Groups") + 
  labs(
    title = "Substance P levels across culture groups",
    subtitle = "N = 3-7, 2 Plates")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = subp1, aes(x = Groups.f, y = sub.p),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot.subp3

```

```{r}
# Subset Values
subp.p1 <- subset(subp0, plate == 1)

m.count.plate1 <- aggregate(data = subp.p1,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
m.count.plate1

# Plot Values
plot.pre.subp0 <- ggplot(subp.p1, aes(x=Groups.f, y=sub.p, fill=Groups.f)) +
  ylab("Substance P [pg/mL]") + 
  xlab("Groups") + 
  labs(
    title = "First quantification plate indicated trends with several zeros present",
    subtitle = "N = 3-7 plate 1")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = subp.p1, aes(x = Groups.f, 
                                y = sub.p),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot.pre.subp0

plot.pre.subp1 <- ggplot(subp.p1, aes(x=Groups.f, y=sub.p)) +
  ylab("Substance P [pg/mL]") + 
  xlab("Groups") + 
  labs(
    title = "Select Samples Identified for Re-Anlaysis",
    subtitle = "N = 3-7 plate 1",
    fill= "Sample")+
  scale_fill_discrete(label=c("No change","Rerun")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red") + 
  geom_dotplot(data = subp.p1, aes(x = Groups.f, 
                                y = sub.p,
                                fill = qual.rem),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1)

plot.pre.subp1

subp.p2 <- subset(subp0, rerun == 1)

plot.pre.subp3 <- ggplot(subp.p2, aes(x=name, y=sub.p, fill=qual.rem)) +
  ylab("Substance P [pg/mL]") + 
  xlab("Sample") + 
  labs(
    title = "Rerun extreme samples were normalized upon reanalysis",
    subtitle = "Negative control = Rat 19 PFF",
    fill= "Measurement")+
  scale_fill_discrete(label=c("Rerun","Original")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "red")

plot.pre.subp3
```



```{r}
# Define Values
subp$subp.4pl <- subp$X4pl.sub.p
subp$subp.4pl[is.na(subp$subp.4pl)] <- 0

# Transform Assessment via histogram
hist(subp$subp.4pl)
hist(subp$t.4pl.subp <- subp$X4pl.sub.p^(1/8))
skewness(subp$t.4pl.subp, na.rm = TRUE)

# Model Development
m.subp.att <- lmer(subp$X4pl.sub.p^(1/8) ~ Groups.f + (1|rat.no), data = subp)
m.subp.att2 <- lmer(subp$X4pl.sub.p^(1/8) ~ Groups.f + (1|rat.no) + (1|plate), data = subp)
m.subp <- lm(subp$X4pl.sub.p^(1/8) ~ Groups.f, data = subp)

# Indentify optimal model
anova(m.subp.att,m.subp,m.subp.att2)

# Statistical summary
em.subp.att <- emmeans(m.subp.att, ~Groups.f)
pwpm(em.subp.att)

em.subp <- emmeans(m.subp, ~Groups.f)
pwpm(em.subp.att)
```